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Abstract 

Using Monte Carlo simulations and mean field calculations for a cell model of water we find 
a dynamic crossover in the orientational correlation time r from non-Arrhenius behavior at high 
temperatures to Arrhenius behavior at low temperatures. This dynamic crossover is independent 
of whether water at very low temperature is charaterized by a "liquid-liquid critical point" or by 
the "singularity free" scenario. We relate r to fluctuations of hydrogen bond network and show 
that the crossover found for r for both scenarios is a consequence of the sharp change in the average 
number of hydrogen bonds at the temperature of the specific heat maximum. We find that the 
effect of pressure on the dynamics is strikingly different in the two scenarios, offering a means to 
distinguish between them. 
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Two different scenarios are commonly used to interpret the anomalies of water 



• The liquid-liquid critical point (LLCP) scenario hypothesizes that supercooled water 
has a liquid-liquid phase transition line that separates a low-density liquid (LDL) at 
low temperature T and low pressure P and a high-density liquid (HDL) at high T and 
P and terminates at a critical point C j^. From C emanates the Widom line Tvi/(P), 
the line of maximum correlation length in the (T, P) plane. Response functions, such 
as the isobaric heat capacity Cp, the coefficient of thermal expansion ap, and the 
isothermal compressibility Kt, have maxima along lines that converge toward Tiy{P) 
upon approaching C [Figs. [1] and [2t^a)]. 

• The singularity-free (SF) scenario hypothesizes the presence of a line of temperatures 
of maximum density Tmd(-P) with negative slope in the (T, P) plane. As a consequence, 
Kt and \ap\ have maxima that increase upon increasing P, as shown using a cell 
model of water. The maxima in Cp do not increase with P, suggesting that there is 
no singularity [4] [Fig. [2](b)]. 

Above the homogeneous nucleation line Th{P) where data are available, the two scenarios 
predict roughly the same equilibrium phase diagram. Here we show that dynamic measure- 
ments should reveal a striking difference between the two scenarios. Specifically, the low-T 
dynamics depends on local structural changes, quantified by the variation of the number 
of hydrogen bonds, that are affected by pressure differently for each scenario. We find this 
result by studying — using Monte Carlo (MC) simulations and mean field calculations — a 
cell model which has the property that by tuning a parameter its predictions conform to 
those of either the LLCP or the SF scenario. This cell model is based on the experimental 
observations that on decreasing P at constant T, or on decreasing T at constant P, (i) 
water displays an increasing local tetrahedrality |5|, (ii) the volume per molecule increases 
at sufficiently low P or T, and (iii) the 0-0-0 angular correlation increases 6j. 

The entire system is divided into cells i G [1, . . . , A^], each containing a molecule with 
volume V = V/N, where V > Nvh_c is the total volume of the system, and fhc is the hard- 
core volume of one molecule. The cell volume f is a continuous variable that gives, in d 
dimensions, the mean distance r = v^^'^ between molecules. The van der Waals interaction 
is represented by a potential with attractive energy e > between nearest-neighbor (n.n.) 
molecules and a hard-core repulsion at i?hc = ^hc*^- 



For a regular square lattice, each molecule i has four bond indices aij G [1, . . . , g], corre- 
sponding to the four n.n. cells j, giving rise to q'^ different molecular orientations. Bonding 
and intramolecular (IM) interactions are accounted for by the two Hamiltonian terms 

7^B = -J$^5..,.,., (1) 

{id) 

where the sum is over n.n. cells, < J < e is the bond energy, 6a,b = 1 if a = 6 and 6a,b = 
otherwise, and 

T-tm = - Ja ^<T,fc<7«, (2) 

where ^^-^ denotes the sum over the IM bond indices {k, I) of the molecule i and Ja > 
is the IM interaction energy with J„ < J, which models the angular correlation between 
the bonds on the same molecule. The total energy of the system is the sum of the van der 
Waals interaction and Eqs. ([T]) and ([2]). 

At constant P, the density of water decreases for T < Tmd(-P) which the model takes 
into account by increasing the total volume by an amount vb > for each bond formed. 
Hence the total molar volume v of the system is 

V = ffree + ^PbVb, (3) 

where vr-ree is a variable for the molar volume without taking into account the bonds, pb 
Nb/{2N) is the fraction of bonds formed and Nb is the number of bonds [4 



We perform simulations in the NPT ensemble j?! for g = 6, vb/v^c = 0.5, J/e = 0.5, and 
for two different values of Jo-/e: (i) Jo-/e = 0.05, which gives rise to a phase diagram with a 
LLCP [Fig. [^a)], and (ii) = 0, which leads to the SF scenario We study two square 
lattices with 900 and 3600 cells, and find no appreciable size effects. We collect statistics 
over 10^ MC steps after equilibrating the system for all P and T. 

For Jcr/e = 0.05, \ap\ for P < Pqi displays a maximum, a™^^ [Fig. [Hb)]. As P increases, 
^max increases and shifts to lower T, converging toward T\y{P) [Fig. [T](a)]. We find that the 
number of bonds, Nb-, increases on decreasing T, and at constant T decreases for increasing 
P, and is almost constant at Tw(-P) Isl. This is consistent with trends seen both in experi- 
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ments [5[ and in simulations [9l|, suggesting that for T > TwiP) the liquid is less structured 
and more HDL-like, while for T < Tw{P) it is more structured and more LDL-like. 

We find that \dpB/dT\ shows a clear maximum for all P < Pc which shifts to lower 
T upon increasing P [Fig. [11(c)] . Remarkably, we also find that the locus of IclpB/dT]'^^'^ 



coincides with the Widom hne Tw{P) [Fig. [T](a)] and that the value of \dpB/dT\^^^ increases 
on approaching Pc- This is the same quahtative behavior as \ap{T)\^^^ and Cp(T)™^^, 
which are used to locate T]y(^P) [Figs. [T]^b) and[2](a)]. The relation of \dpB/dT\ with the 
fluctuations is revealed by its proportionality to |a;p(T)| and to the fluctuation of the number 
of bonds 



2NkBT^ 



dp 



B 



dT 



(4) 



J~PVB 

where ks is the Boltzmann constant. 

For Jo- = (SF scenario) we observe no difference for the behavior of Nb and \dp-B/dT\. 
We further verify the prediction of the SF scenario ^ that Cp'^^ remains a constant upon 
increasing P [Fig. [2](b)]. 

Next, we study how this different behavior affects the dynamics. Previous simulations 
[lo| found a crossover from non-Arrhenius to Arrhenius dynamics for the diffusion constant 
of models that display a LLCP, and showed the temperature of this crossover coincided with 
T]v{P). We calculate, for both scenarios, the relaxation time r of Si = Ylj'^ijl^i "which 
quantifles the degree of total bond ordering for site i. Speciflcally, we identify r as the time 
for the spin autocorrelation function C„^{t) = {Si(t)Si{0)) / {Sf{0)) to decay to the value 
1/e. 

For both scenarios we flnd a dynamic crossover (Fig. [3]). At high T, we flt r with the 
Vogel-Fulcher-Tamman (VFT) function 

T-To 

where Tq^^, Ti, and To are three fltting parameters. We flnd that r has an Arrhenius T 
dependence at low T, t = tq exp[i?A/^_BT], where tq is the relaxation time in the high- 
T limit, and Ea is a T-independent activation energy. We flnd that for Jo-/e = 0.05 the 
crossover occurs at T]v{P) for P < Pq' [Fig. [3]^a)], and that for J^- = the crossover is at 
T(Cp'^^), the temperature of C^^^ [Fig. [3]^b)]. We note that for both scenarios the crossover 
is isochronic, i.e. the value of the crossover time tq is approximately independent of pressure. 
We flnd Tc ~ 103/2 MC steps. 

We next calculate the Arrhenius activation energy Ea{P) from the low-T slope of logr 
vs. 1/T [Fig. 111(a)]. We extrapolate the temperature Ta(P) at which r reaches a flxed 



VFT _ VFT 

T = Tq exp 



(5) 



macroscopic time ta > tq. We choose ta = 10 MC steps > 100 sec 



Ul [Fig.HKb)]. We 



flnd that Ea{P) and Ta{P) decrease upon increasing P in both scenarios, providing no 



distinction between the two interpretations. Instead, we find a dramatic difference in the 
P dependence of the quantity E^/ {kBTjC) in the two scenarios, increasing for the LLCP 
scenario and approximately constant for the SF scenario [Fig. ID^c)]. 

We can better understand our findings by developing an expression for r in terms of 
thermodynamic quantities, which will then allow us to explicitly calculate Ep^/ [ksTp^ for 
both scenarios. For any activated process, in which the relaxation from an initial state to a 
final state passes through an excited transition state, ln(r/ro) = A(f/ + PV — TS)/{kBT), 
where A(U + PV — TS) is the difference in free energy between the transition state and the 
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initial state. Consistent with results from simulations and experiments p^, 113|], we propose 
that at low T the mechanism to relax from a less structured state (lower tetrahedral order) to 
a more structured state (higher tetrahedral order) corresponds to the breaking of a bond and 
the simultaneous molecular reorientation for the formation of a new bond. The transition 
state is represented by the molecule with a broken bond and more tetrahedral IM order. 
Hence, 

A{U + PV- TS) = JpB - JaPiM - PvB - TAS, (6) 

where pb and piu, the probability of a satisfied IM interaction, can be directly calculated. 
To estimate AS, the increase of entropy due to the breaking of a bond, we use the mean 
field expression AS = fcs[ln(2A^pB) — ln(l + 2N{1 — Pb))]pb, where pB is the average value 
of pb above and below Tw{P)- 

We next test that the expression of ln(r/ro), in terms of AS and Eq.([n]), 

T JpB - JaPlM - Pvb - 1 ^NpB 

To KbT 1 + 2A^(1 - Pb) 

describes the simulations well, with minor corrections at high T. Here tq = to{P) is a free 
fitting parameter equal to the relaxation time for T oo. From Eq.(I71) we find that the 
ratio Ex/ {ksTx) calculated at low T increases with P for J„/e = 0.05, while it is constant 
for Jo- = 0, as from our simulations [Fig. Hl^d)]. 

In summary, we have seen that both the LLCP and SF scenarios exhibit a dynamic 
crossover at a temperature close to T(Cp'*^), which decreases for increasing P. We interpret 
the dynamic crossover as a consequence of a local breaking and reorientation of the bonds 
for the formation of new and more tetrahedrally oriented bonds. Above T{Cp^^), when T 
decreases, the number of hydrogen bonds increases, giving rise to an increasing activation 
energy E'a and to a non-Arrhenius dynamics. As T decreases, entropy must decrease. A 



major contributor to entropy is the orientational disorder, that is a function of pb, as de- 
scribed by the mean field expression for AS. We find that, as T decreases, pb — hence 
the orientational order — increases. We find that the rate of increase has a maximum at 
T(Cp'^^), and as T continues to decrease this rate drops rapidly to zero — meaning that 
for T < T{Cp^^), the local orientational order rapidly becomes temperature-independent 
and the activation energy Ea also becomes approximately temperature-independent, for the 
Eq-©. Corresponding to this fact the dynamics becomes approximately Arrhenius. 

We find that the crossover is approximately isochronic (independent of the pressure) 
consistent with our calculations of an almost constant number of bonds at T{Cp ''^^). In both 
scenarios, and Ta decrease upon increasing P, but the P dependence of the quantity 
Ea/ {^bTa) has a dramatically different behavior in the two scenarios. For the LLCP scenario 
it increases as P — > Pc"? while it is approximately constant in the SF scenario. We interpret 
this difference as a consequence of the larger increase of the rate of change of pb in the 
LLCP scenario, where pb diverges at finite T^/, compared to the SF scenario, where pb 
can possibly diverge only at T = 0. Since experiments can detect local changes of water 



structure from HDL-like to LDL-like, (e. g. IJ]), it is possible that our prediction on the 



dynamic consequences of this local change may be experimentally testable. 
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FIG. 1: (color online) (a) The relevant part of the liquid phase diagram for the LLCP scenario. 
C is the LDL-HDL critical point at the end of the first-order phase transition line (thick line) Q], 
Tvi/(P) (thin line) is the Widom line, which we take to be the average between T([ap|™'^^) (Q) a-nd 
T(Cp'*^) (□) from panel (b), and Fig. [^Ja), respectively. Upper and lower dashed line are quadratic 
fits of |ap|™'^^ and Cp*^^, respectively, merging at the C . T(|ap|™'^^) and T(Cp'^^) overlap within 
the error bars, (b) \ap\ as a function of T for different P. (c) [dpe/dTl, the temperature derivative 
of pB, as a function of T for different P. In panel (a), r(|(ipB/c?P|™^^) (^) overlaps Tt^{P\ 



8 




FIG. 2: (color online) Temperature dependence of the specific heat Cp for both the LLCP and 
the SF scenarios, for seven values of P. (a) For the LLCP scenario, Cp has a maximum, the size 
of which increases with increasing pressure and diverges as P ^ Pq' ■ (b) For the SF scenario, Cp 
also has a maximum, but its size does not increase with increasing pressure, consistent with the 
findings of the mean field calculations of Ref . jj] . 




e/k,T 



FIG. 3: Dynamic crossover — large hatched circles of a radius approximately equal to the error 
bar — in the orientational relaxation time r for a range of different pressures, (a) The LLCP 
scenario, with crossover temperature at Tw[P). (b) The SF scenario, with crossover temperature 
at r(Cp^^). Solid and dashed lines represent Arrhenius and VFT fits, respectively. Notice that 
the dynamic crossover occurs at approximately the same value of r for all seven values of pressure 
studied. 
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FIG. 4: Effect of pressure on the activation energy Ex- (a) Demonstration that Ep^ decreases 
linearly for increasing P for both the LLCP and the SF scenarios. The lines are linear fits to 
the simulation results (symbols), (b) Ta, defined such that t{Ta_) = W^'^ MC steps > 100 sec 



11[, decreases linearly with P for both scenarios, (c) P dependence of the quantity E'a/C^bTa) is 
different in the two scenarios. In the LLCP scenario, E^/ {kBTx) increases with increasing P, and 
it is approximately constant in the SF scenario. The lines are guides to the eyes, (d) Demonstration 
that the same behavior is found using the mean field approximation. In all the panels, where not 
shown, the error bars are smaller than the symbol sizes. 
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